import pulp
import pandas as pd
from matplotlib import font_manager
import matplotlib.pyplot as plt

# 配置中文字体显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 读取用户上传的数据
file_path = '附件1.xlsx'
df_land = pd.read_excel(file_path, sheet_name='乡村的现有耕地')

file_path_2 = '附件2.xlsx'  # 确保路径指向正确的文件
df_crops = pd.read_excel(file_path_2, sheet_name='2023年统计的相关数据')

# 定义豆类作物的编号集合
bean_crops = {1, 2, 3, 4, 5, 17, 18, 19}

# 整理地块类型并分类
A = {row['地块名称']: row['地块面积/亩'] for _, row in df_land.iterrows()}

# 单季地块：平旱地、梯田、山坡地
plots_single_season = df_land[df_land['地块类型'].isin(['平旱地', '梯田', '山坡地'])]['地块名称'].tolist()

# 双季地块：水浇地、普通大棚、智慧大棚
plots_double_season = df_land[df_land['地块类型'].isin(['水浇地', '普通大棚', '智慧大棚'])]['地块名称'].tolist()
#将双季节再进行细分一下
irrigated_land = df_land[df_land['地块类型'] == '水浇地']['地块名称'].tolist()
greenhouse_land = df_land[df_land['地块类型'] == '普通大棚']['地块名称'].tolist()
smart_greenhouse_land = df_land[df_land['地块类型'] == '智慧大棚']['地块名称'].tolist()
# 创建问题实例
model = pulp.LpProblem("Crop_Planting_Optimization", pulp.LpMaximize)

# 定义作物类型
crops_single_season = list(range(1, 16))  # 作物1-15，适用于单季地块
crops_double_season_irrigated = list(range(16, 38))  # 作物16-37，适用于水浇地
crops_double_season_greenhouse = list(range(23, 42))  # 作物23-41，适用于普通大棚
crops_double_season_smart = list(range(27, 35))  # 作物27-34，适用于智慧大棚
df_crops['销售单价(元/斤)'] = df_crops['销售单价/(元/斤)'].apply(
    lambda x: (float(x.split('-')[0]) + float(x.split('-')[1])) / 2 if isinstance(x, str) else x
)

# 生成销售单价、亩产量和种植成本的字典
P = {row['作物编号']: row['销售单价(元/斤)'] for _, row in df_crops.iterrows()}  # 销售单价字典
Y = {row['作物编号']: row['亩产量/斤'] for _, row in df_crops.iterrows()}  # 亩产量字典
C = {row['作物编号']: row['种植成本/(元/亩)'] for _, row in df_crops.iterrows()}  # 种植成本字典

# 定义决策变量 x[i][k][t] 是地块 i 在年份 t 的季节 j 是否种植作物 k（取值为0, 0.5或1）、
years = list(range(2024, 2031))  # 从2024到2030
x = pulp.LpVariable.dicts("x", (plots_single_season + plots_double_season,
                                crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart,
                                years, [1, 2]), 0, 1, cat='Continuous')
# 定义辅助二进制变量 y1 和 y2
y1 = pulp.LpVariable.dicts("y1", (plots_single_season + plots_double_season,
                                  crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart,
                                  years, [1, 2]), 0, 1, cat='Binary')

y2 = pulp.LpVariable.dicts("y2", (plots_single_season + plots_double_season,
                                  crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart,
                                  years, [1, 2]), 0, 1, cat='Binary')

# 目标函数：最大化所有地块的净收益
# 单季地块目标函数
# 单季地块目标函数
Z_single = pulp.lpSum(
    (P[k] * Y[k] * A[i] - C[k] * A[i]) * x[i][k][t][1] 
    for t in years 
    for i in plots_single_season 
    for k in crops_single_season
)

# 双季地块目标函数
Z_irrigated = pulp.lpSum(
    (P[k] * Y[k] * A[i] - C[k] * A[i]) * (x[i][k][t][1] + x[i][k][t][2])
    for t in years
    for i in irrigated_land
    for k in crops_double_season_irrigated
)

Z_greenhouse = pulp.lpSum(
    (P[k] * Y[k] * A[i] - C[k] * A[i]) * (x[i][k][t][1] + x[i][k][t][2])
    for t in years
    for i in greenhouse_land
    for k in crops_double_season_greenhouse
)

Z_smart_greenhouse = pulp.lpSum(
    (P[k] * Y[k] * A[i] - C[k] * A[i]) * (x[i][k][t][1] + x[i][k][t][2])
    for t in years
    for i in smart_greenhouse_land
    for k in crops_double_season_smart
)

# 综合目标函数
model += Z_single + Z_irrigated + Z_greenhouse + Z_smart_greenhouse
#x变量约束条件
for i in plots_single_season + irrigated_land + greenhouse_land + smart_greenhouse_land:
    for k in crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart:
        for t in years:
            for j in [1, 2]:
                # 约束 x 的值等于 0.5 * y1 + y2
                model += x[i][k][t][j] == 0.5 * y1[i][k][t][j] + y2[i][k][t][j]
                # 确保 y1 和 y2 不会同时为 1
                model += y1[i][k][t][j] + y2[i][k][t][j] <= 1
# 约束条件1：为每种地块类型添加作物选择约束
for t in years:
    # 单季地块（平旱地、梯田、山坡地）限制作物选择
    for i in plots_single_season:
        # 在每年只能从作物 {1, 2, ..., 15} 中选择
        model += pulp.lpSum(x[i][k][t][1] for k in crops_single_season) == 1
    # 双季地块
    # 水浇地
    for i in irrigated_land:
        # 第一季可以在 {16, 17, ..., 34} 中选择
        model += pulp.lpSum(x[i][k][t][1] for k in range(16, 35)) == 1
        # 第二季只能在 {16, 35, 36, 37} 中选择
        model += pulp.lpSum(x[i][k][t][2] for k in [16, 35, 36, 37]) == 1

    # 普通大棚
    for i in greenhouse_land:
        # 第一季可以在 {17, 18, ..., 34} 中选择
        model += pulp.lpSum(x[i][k][t][1] for k in range(17, 35)) == 1
        # 第二季只能在 {38, 39, 40, 41} 中选择
        model += pulp.lpSum(x[i][k][t][2] for k in range(38, 42)) == 1

    # 智慧大棚
    for i in smart_greenhouse_land:
        # 第一季和第二季都只能在 {17, 18, ..., 34} 中选择
        model += pulp.lpSum(x[i][k][t][1] for k in range(17, 35)) == 1
        model += pulp.lpSum(x[i][k][t][2] for k in range(17, 35)) == 1

# 约束条件2：在每块地上每季的作物种植比例之和为1
for t in years:
    for i in plots_single_season:
        model += pulp.lpSum(x[i][k][t][1] for k in crops_single_season) == 1

    for i in plots_double_season:
        for j in [1, 2]:
            model += pulp.lpSum(x[i][k][t][j] for k in crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart) == 1

# 约束条件3：轮作约束
for t in years[:-1]:
    for i in plots_single_season:
        for k in crops_single_season:
            model += x[i][k][t][1] + x[i][k][t+1][1] <= 1

    for i in plots_double_season:
        for k in crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart:
            model += x[i][k][t][1] + x[i][k][t][2] <= 1
            model += x[i][k][t][2] + x[i][k][t+1][1] <= 1

#越苏条件3：三年内需要存在豆类

# 定义约束：每块地在每连续三年内至少种植一次豆类作物
for i in plots_single_season:
    for t_start in range(2024, 2028):  # 每个三年窗口的起始年份（最多到2027）
        model += pulp.lpSum(x[i][k][t][1] for k in bean_crops for t in range(t_start, t_start + 3)) >= 1

for i in irrigated_land + greenhouse_land + smart_greenhouse_land:
    for t_start in range(2024, 2028):  # 每个三年窗口的起始年份（最多到2027）
        model += pulp.lpSum(x[i][k][t][1] + x[i][k][t][2] for k in bean_crops for t in range(t_start, t_start + 3)) >= 1

# 求解模型
model.solve()

# 创建一个空列表来存储结果
results = []

# 遍历所有变量，并将非零决策变量的值保存下来
for i in plots_single_season + plots_double_season:
    for k in crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart:
        for t in years:
            for j in [1, 2]:
                if x[i][k][t][j].varValue is not None and x[i][k][t][j].varValue > 0:
                    results.append([i, k, t, j, x[i][k][t][j].varValue])

# 使用 pandas 将结果转换为数据框
df_results = pd.DataFrame(results, columns=["地块", "作物", "年份", "季节", "种植比例"])

# 显示结果
print(df_results)
output_path = "optimization_results4.csv"
df_results.to_csv(output_path, index=False)
print("保存成功！")



import matplotlib.pyplot as plt
import seaborn as sns

# 设置绘图的风格
sns.set(style="whitegrid")
# 将 matplotlib 的默认字体修改为黑体，字号修改为12
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定英文字体为Times New Roman
plt.rcParams['axes.unicode_minus'] = False    # 解决负数坐标轴显示问题
plt.rcParams['font.size'] = 15                # 指定字号

# 数据准备，获取不同年份、地块和作物种植比例的汇总
df_yearly = df_results.groupby(["年份", "作物"]).agg({"种植比例": "sum"}).reset_index()

# 图表1：作物种植面积随时间的变化图
plt.figure(figsize=(10, 6))
for crop in df_yearly["作物"].unique():
    crop_data = df_yearly[df_yearly["作物"] == crop]
    plt.plot(crop_data["年份"], crop_data["种植比例"], label=f"作物 {crop}", marker='o')

plt.title("各种作物种植面积随时间的变化", fontsize=14)
plt.xlabel("年份", fontsize=12)
plt.ylabel("种植比例", fontsize=12)
plt.legend(loc="upper right", title="作物名称")
plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt
import pulp
import pandas as pd

# 记录迭代过程中的目标函数值
objective_values = []
iterations = []

def capture_objective_value():
    # 计算当前模型的目标函数值
    current_value = pulp.value(model.objective)
    objective_values.append(current_value)
    iterations.append(len(objective_values))  # 记录当前是第几次迭代

# 修改模型求解过程，设置一个回调函数来捕捉每次迭代的目标函数值
def solve_with_callback(model):
    # 在求解之前增加捕获函数
    model.solve(pulp.PULP_CBC_CMD(msg=True, warmStart=True))
    capture_objective_value()

# 模型创建部分保持不变
# ... (省略之前的模型定义和约束部分)

# 求解模型并记录每次迭代
solve_with_callback(model)

# 创建一个空列表来存储结果
results = []

# 遍历所有变量，并将非零决策变量的值保存下来
for i in plots_single_season + plots_double_season:
    for k in crops_single_season + crops_double_season_irrigated + crops_double_season_greenhouse + crops_double_season_smart:
        for t in years:
            for j in [1, 2]:
                if x[i][k][t][j].varValue is not None and x[i][k][t][j].varValue > 0:
                    results.append([i, k, t, j, x[i][k][t][j].varValue])

# 使用 pandas 将结果转换为数据框
df_results = pd.DataFrame(results, columns=["地块", "作物", "年份", "季节", "种植比例"])

# 保存结果
output_path = "optimization_results4.csv"
df_results.to_csv(output_path, index=False)
print("保存成功！")

# 图表1：模型目标函数随迭代次数的变化图
plt.figure(figsize=(10, 6))
plt.plot(iterations, objective_values, marker='o', color='blue', linewidth=2)
plt.title("模型目标函数随迭代次数的变化", fontsize=14)
plt.xlabel("迭代次数", fontsize=12)
plt.ylabel("目标函数值 (元)", fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()

# 图表2：决策变量随迭代次数的变化图
# 假设我们选择某个关键决策变量，比如某个地块的种植面积变化
# 这里使用决策变量变化的实际值
key_decision_variable = [x[i][1][2024][1].varValue for i in plots_single_season]  # 使用实际数据

plt.figure(figsize=(10, 6))
plt.plot(iterations, key_decision_variable, marker='o', color='red', linewidth=2)
plt.title("决策变量随迭代次数的变化", fontsize=14)
plt.xlabel("迭代次数", fontsize=12)
plt.ylabel("种植面积比例", fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()
